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ABSTRACT 



The general equations of surface ship motion are 
developed and standardized for simulation in digital 
computer. Digital simulations of the dynamics of the 
surface ship in three degrees of freedom are done 
with and without non-linear and cross-coupling terms. 
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I . EQUATIONS OF SURFACE SHIP MOTION 



A moving ship is a body with six degrees of freedom. 
These degrees of freedom are generally chosen as follows: 

a. Linear displacements along the three axes 
through the center of gravity . 

a.l Surge — along X axes 
a. 2 Sway — along Y axes 
a. 3 Heave — along Z axes 

b. Rotations around the three axes through the 
center of gravity. 



b.l 


Roll — 


around 


X axes 


b . 2 


Pitch — 


around 


Y axes 


b . 3 


Yaw 


around 


Z axes 


Further 


reduction 


in the 


complex nature of the equa 



tions can be brought about by choosing an orthogonal 
axis system parallel to the principal axes of inertia so 
as to eliminate products of inertia in the motion equa- 
tions. For practically all ocean vehicles, with extreme- 
ly few exceptions, a longitudinal axis (X axis) in the 
centerline plane, a downward (toward keel) axis (Z axis) 
perpendicular to the X axis in the centerline plane, 
and a transverse axis (Y axis) perpendicular to the 
centerline plane satisfies this requirement. 
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SUK&Z + 



Figure 1. Surface ship in six degrees of freedom. 



For the exceptional vehicle which has a very peculiar 
and significantly large asymmetrical mass distribution, 
it is necessary to include the products of inertia. 

The X, Y, Z axes form an orthogonal right hand system 
of axes fixed in the vehicle. The axes and the associated 
components of the pertinent physical quantities are 
defined below: 

The longitudinal X axis (in the plane of symmetry) is 
positive in the forward direction, usually parallel to 
the keel or calm water line. If the upper and lower halves 
of the body are symmetrical, then the axis is the inter- 
section of the two planes of symmetry. 

The Y axis is the transverse axis perpendicular to the 
plane of symmetry and positive to the starboard. 

The Z axis or downward axis in the plane of symmetry 
(X,Z) is perpendicular to the X axis and positive downward 
towards the keel. 
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AAA 

i,j,k unit vectors along the X,Y,Z axis respectively. 

R x , y , z vector distance of a point from the origin 0, 
and the corresponding components along the X,Y and Z axes. 

R = ix G + jy G + kz G 

U u,v,w velocity or the origin 0 (on the body) and 
the corresponding components along the X,Y and Z axis. 

U = iu + jv + kw 

ti p,q,r angular velocity of the body about the origin 
and the corresponding components about the axes. 

A A A 

H = ip + jn + kr 

The moments of inertia of the body about the X,Y and 
Z axes respectively 1^, 1^ , 1^. 

F X , Y , Z , force acting on the body and the correspond- 
ing components along the axes. 

^ A A A 

F = iX + jY + kZ 

m K,M,N Moments acting about the axes. 

^ A A A 

m = iK + jM + kN 

Newton's law of motion for a rigid body can be written 
as two equations, one a force equation and the second a 
moment equation provided an origin is taken at the center 
of gravity and the axis system is fixed in space. The 
equations are 

d > d 

F = dt ( mornenturn ) * ( m U ) 

d > d 

m = (angular momentum) = g^-(I^) 
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where the subscript G refers to an origin at the center 
of gravity and m is the mass of the body. For a mass 
essentially constant in time 



F 




G 



) 



For an origin not at the center of gravity of the 
body and in a system of axes fixed in and moving with 
the vehicle. 

U_ = U + ft x R_ 
a. o 

where U is the velocity of the origin in space. However, 

cl 

since the origin is on the surface of the earth and the 
earth rotates, then 



— >■ -V — ^ 

U = U + fl x R, 
a e b 

where U is the geographical velocity of the body, is 
the angular velocity of the earth, and is the radius 
vector from earth's center to the vehicle. The force 
equation becomes: 



cl . ->- ->* 

F = m g^-(U + !2 g x x R^) 

F = m gjr(U + ft x R^J+mt^ x R^+Q^x R^] 



= m 



ir-(U + & x R„)+m[^ x U x ^ x R, ] 
flt G e e e b 



since 



ti = 0 and R^ = U + 



*b 



the term mft xU is the coriolis force and mS^ x x R, 
e e e b 

is the centripetal acceleration due to rotation of the 
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earth. These two terras are negligibly small when com- 
pared with the other forces, then 



F = m x R q ) 

Finding the derivatives of unit vectors (change in 
direction) 
where 



A 

di = 


-kd 0 


A 

di 


-=>- 

< 

II 




A 

di = 


= 0 


/s 

dj = 


0 


A 

dj 


= -id t 




A 

dj = 


= kd|) 


/s 

dk = 


id 9 


A 

dk 


= 0 




A 

dk = 


= -jd* 


the 


contributions 












/N 

di 

dt 


/s 

i • 0 + 


- £ 
3 dt K 


‘d0 

dt 








A 

dj _ 
dt 


+ 

dt 


A 

j • 0 + 


k d ^ 

K dt 








/\ 

dk 


/ d0 


Cd(j) 


A 








dt 


1 dt 


3 dt + 


k • 0 







. A A A 

ft = ip + jq + kr 
_ d(J> d0 



p dt ' q “ dt ' 



r = §1 
dt 



di 

dt 



= a 



X 1 = 



1 

p 



3 

q 

0 



k 

r 

0 



"> , t 
U = iu + jv + kw 

R G = iX G + ^ Y G + kZ G 



dU 



A A J ' A T • A 

= id + u(j2_ + jv + v^J- + kw + w 



dk 



dt ' ~ VJ dt ' JV ' v dt 1 1 "dt 

- ii + U(jS| - k |f-) + jV+ V(-i|| + k||> 



+ ^ + .«i|| - in, 



” i ” ^ A ^ A A A 

= iu+ u(jr- kq)+ jv+ v (~ir+ kp)+kw + w(iq - jp) 
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3 



x r g = 



yN /\ S\ 

i j k 

p q r 



X G Y G Z G 



= i (qZ - rY )- j (pZ - rX )+ k (pY -qX ) 



d -(^XR G ) = i(Z G q-Y G r)+ j (Z G qr-Y G r 2 ) + k(Y G qr-Z G q 2 ) 



dt 



A A A A A 

+ j (X G r - Z G P)+ k (X Q rp - Z Q p )+ i(Z G rp - X Q r ) 

+ k(Y G p - X Q q)+ i (Y Q pq - X Q q 2 )+ j (X Q pq - Y Q p 2 ) 

from F = iX + jY + kz and F = m[^- + ^(£2 x R Q ) 

/v. /\ 

Rewriting and grouping all i terms equal to X, all j 

A 

terms equal to Y and all k terms equal to Z yield 



X = m (u-vr+wg+Z^q + rpZ Q -r 2 X G + pqY G - q 2 X G ) 



G 

C G 

C G 

C G 



2 2 

= m (u+wg-vr-X^, (r +q ) + Y G (pq-r) + Z G (rp+q) 

Y = ra (v+ur-wp+X^, (pq+r ) - Y G (p 2 +r 2 ) + Z G (pq-r) 

Z = m (w+vp-uq+X^, (rp-q) + Y G (qr+p) - Z Q (p 2 + q 2 ) 



( 1 ) 



From 



m G = ^-(angular momentum) 
d 



(il v P + j I v <3 + kI «7 r ) 



dt ' X 



G indicates an origin at the center of gravity 
I X = I X “ m (Y G + Z G ) 

J Y g ~ I Y ■ 10 <Z G + X G 1 
“ (X G 2 + Y G 2) 

M = M g + Rq X F or M g = M - R g X F 
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After manipulating and using the results for the 
derivatives of unit vectors (same as above) expressions 
for K, H and N are obtained. 

K = I x p +(I Z - I y ) qr+m[Y Q (w+ pv - qu) - Z Q (v+ru-pw) ] 

M = I y q +( I X - I z )rp+m[Z G (u+ qw - rv) - X Q (w+pv-qu) ) (2) 

N = I z r +(I y - I x )pq+m[X G (v+ ru - pw) - Y G (u+qw-rv)) 

The terms (qw-rv) , (ru-pw) and (pv-qu) are gyroscopic 

effects . 

The relationship for forces and moments can be ex- 
pressed 

Force) 

) = f (properties of body, properties of 
Moment) motion, properties of fluid) 



For a particular ship, in a given fluid with no 
excitation force - so 



F ) 

^ ) = f (properties of motion) 

M ) 

= f(X o ,Y o ,Z o ,<j),0,^,u,v,w,p,q,r,u,v,w, 
p,q,r,<$,6 ' ,6") 



The Taylor series which has the following form may 
now be applied to linearize the equations about an 
operating point X Q 



f(X) = f(X Q ) + 



(x - x o>s s<v+<x-x o ) 2 r-2 

dX 



f(X )+ 

o 



Apply this to f(X,Y,Z) 
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For the case of f (properties of motion) let (X-X ) 
= AX q , (Y-Y o ) = AY q and (Z-Z Q ) = AZ q 

terms second order and higher are neglected for small 
perturbations . 

From X equation the linear terms are obtained: 



X = f( ) + AX (|| )+ AY q (|| )+ AZ q (|| )+ . 

o o o 



A / 6 f \ , 

•* Av( ^7 )+ * 



The defining relations are: 



/ <5 f \ = 

^3"u o ^3"u u=u 



= X 



u 



( 6X ) = X 

'<!Tw w=w w 



Aw = (w-w )=w, w = 0 
o o 



Au = (u-u q ) 



The force equations then become 



X = X 

o 


+ 


X v X 
x o ° 


+ 


0 

0 

X 


+ X Q 0 


+ 


••••Vu 




Y = Y 

o 


+ 


Y X X o 
o 


+ 


Y v Y 
Y o 
o 


+ Y 0 0 


+ 




(3) 


Z = Z 

o 


+ 


Z X X o 
o 


+ 


Y 

Y o 
O 


+ . . . . Z 0 0 


+ 


Z u A u 





A similar derivation can be done for the K,M and N equa- 
tions. The preceding X,Y and Z equations may now be 
equated to the linearized X,Y and Z (equations (1)), e.g., 
for the Y equation without roll, pitch, and the center 
of gravity at X G = 0, Y^= 0 and Z Q = 0 gives the linearized 
equation . 

Y = m(v + ur) 
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then 



Y„ X + Y„ Y + Y . Y u+ Y v+ Y r+Y*u+ Y*v+ Y*r = m(v+ru) 
XoYoiKuvruvr 
o o Y 

Expressions for X,Z and K,M and N can be determined in 
a similar procedure. 

In order to obtain the equations in a non-dimensional 
form some definitions will be given, and applied to the Y 
force equations as an example of the process. 



Froude number 

U 1 .V ' ,w ' 



t' 



X ' , Y ' , Z ' 



U 

/gT 

u,v,w 

/g/T 

t (/g Jl) 

X, Y,Z 

pg ^ 3 



c x’' c y’' c z 



_ X , Y , Z 
1 „2 „ 2 
2 pU £ 



After replacing and adding the effect of waves, the 
Y equation becomes: 

v'+ r'u' = \ U 2, (C v' + C v p' + C r'+ c 6r ) 

V Y p' Y r ' Y 6r 

+ (Y^.p’ + Y^ v’ + Y^.r’ + Y w ') 
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II. DIGITAL COMPUTER SIMULATION 



The six equations of motion after rearranging by plac- 



ing the second order terms to the left and the rest on the 



right side become: 

.. .. .. .. .... • • 

aaA+baB+caC+daD+eaE+f aF = - (a ^a^A+a 2 a 2 A+b^a ^B+b 2 a 2 B 

• • 

+ c^a^C+c 2 a 2 C+d^a^D+d 2 a 2 D 

• • 

+ e^a^E+e 2 a 2 E+f ^a^F+f 2 a 2 F) 
+ IFl 



abA+bbB+cbC+dbD+ebE+fbF = - (a^b^A+a 2 b 2 A+b^b^B+b 2 b 2 B 

* • 

+ c 1 b 1 C+c 2 b 2 C+d 1 b 1 D+d 2 b 2 D 

+ e lV +e 2 b 2 E+£ l b l F+f 2 b 2 F> 

+ IF2 



• • • • •• •• It II ^ • 

acA+bcB+ccC+dcD+ecE+fcF = - (a^c^A+a 2 c 2 A+b^c^B+b 2 c 2 B 

• . 

+ c 1 c 1 C+c 2 c 2 C+d^c 1 D+d 2 c 2 D 

• . 

+ e i c ^E+e 2 c 2 E+f ^c^F+f 2 c 2 F) 
+ IF3 

*. .. •• .. . . .. * • 

adA+bdB+cdC+ddD+edE+f dF = - (a^-jA+a^^+b-^B+b^^ 

♦ * 

+ c^d^C+c 2 d 2 C+d^d^D+d 2 d 2 D 
+ e 1 d 1 E + e 2 d 2 E + f 1 a 1 F + f 2 d 2 F) 
+ IF4 

*• •• .. .. • . .. . . 

aeA+beB+ceC+deD+eeE+feF = - (a^e^A+a 2 e 2 A+b^e^B+b 2 e 2 B 

+ c ]^i c+C 2 e 2 C+d l e l D+d 2 e 2 D 
+ e 1 e 1 E+e 2 e 2 E+f 1 e 1 F+f 2 e 2 F) 

+ IF5 
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/ 




af A+bf B+cf C+dfD+ef E+f fF = - (a 1 f x A+a 2 f 2 A+b 1 f 1 B+b 2 f 2 B 

+ c 1 f 1 C+c 2 f 2 C+d 1 f 1 D+d 2 f 2 D 

+ e l £ l E+e 2 £ 2 E+£ l £ l F+f 2 £ 2 F) 

+ IF6 

•• . ** . * " . * ” . 

where A=u, A=u, B=v, B=v , C=w , C=w, D=p, D=p 

•• . •• . 

E=q, E=q , F=r, F=r , terms IF include all non-linear 



terms such as wave force, wind, rudder deflection etc. 

In the six equations, the right can be set equal to 



I I' I 2 ' 1 6 respectively, thus: 



Z 1 


= 


- (a-^a-^A+a 2 a 2 A+b^a^B+ 


• £ 2 a 2 F> 


+ 


IFl 


J 2 


= 


~ ^ a l^l A+a 2^2 A+ ^ > l^l B+ 


f 2 b 2 F) 


+ 


IF2 


Z 3 


= 


- (a 1 o 1 A+a 2 c 2 A+b 1 c 1 B+ 


f 2 c 2 F) 


4- 


IF3 


*4 


= 


- (a, d, A+a 0 d„A+b., d, B+ 


f 2 d 2 F) 


+ 


IF4 


J 5 


= 


~ ( a i e i A+a 2 e 2 A+B l e l B+ 


f 2 e 2 F ^ 


4- 


IF5 


J 6 




-<a 1 f 1 A+a 2 f 2 A+b 1 f 1 B+ 


f 2 f 2 F) 


4- 


IF6 



the equations then have 



aaA 


4* 


baB 


4* 


caC 


+ 


daD 


+ 


abA 


+ 


bbB 


4* 


cbC 


4* 


dbD 


+ 


acA 


+ 


bcB 


+ 


ccC 


4* 


dcD 


+ 


adA 


+ 


bdB 


4* 


cdC 


4* 


ddD 


4* 


aeA 


+ 


beB 


+ 


ceC 


+ 


deD 


+ 


afA 


+ 


bf B 


4- 


cfC 


+ 


dfD 





the form that follows 

eaE + faF = 1^ 

ebE + fbF = I 

ecE + fcF = I^ 

edE + fdF = I, 

4 

eeE + feF = I c 
o 

ef E 4* f f F = 1^ 

6 



/ 
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expressing in matrix form: 



aa 

ab 

ac 

ad 

ae 

af 



A = 



ba 


ca 


da 


ea 


fa 




A 


bb 


cb 


db 


eb 


fb 




B 


be 


cc 


dc 


ec 


f c 




C 


bd 


cd 


dd 


ed 


fd 




D 


be 


ce 


de 


ee 


fe 




E 


bf 


cf 


df 


ef 


ff 




F 


amer 


' s rule 


to 


solve 


for A, 


B ■ 




T 


1^ 










i 


6 










*1 


ba 


ca 


da 


ea 


fa 


X 2 


bb 


cb 


db 


eb 


fb 


Z 3 


be 


cc 


dc 


ec 


fc 


*4 


bd 


cd 


dd 


ed 


fd 


*5 


be 


ce 


de 


ee 


fe 


*6 


bf 


cf 


df 


ef 


ff 


aa 


ba 


ca 


da 


ea 


fa 


ab 


bb 


cb 


db 


eb 


fb 


ac 


be ' 


cc 


dc 


ec 


f c 


ad 


bd 


cd 


dd 


ed 


fd 


ae 


be 


ce 


de 


ee 


fe 


af 


bf 


cf 


df 


ef 


ff 



F in 



define the denominator determinant _ A and for the 
nominator let cofactor of 1^ = cof. aa 



cofactor of = cof. ab , -■ 
= cof. af equations becomes 



-and coefactor of I, 



A = 



(cof . aal^+cof . abH^+cof . acl^+cof . adl^+cof . ael^+cof . af Ig ) 
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In the same way solve for B, C, 



(cof 


bal 1 


cof 


bbi 2 


cof 


bcl 3 


cof 


bdl 4 


cof 


bel 5 


cof 


bfl 6 ) 


(cof 


cal^ 


cof 


cbl 2 


cof 


A 

ccl 3 


cof 


cdl 4 


cof 


cel,. 

D 


cof 


cf V 


(cof 


dal^ 


cof 


dbl 2 


cof 


A 

dcl 3 


cof 


ddl 4 


cof 


del c 

o 


cof 


dfl 6 ) 


(cof 


eal-^ 


cof 


ebl 2 


cof 


A 

ecI 3 


cof 


edl 4 


cof 


eel- 

D 


cof 


e «6> 


(cof 


fal. 


cof 


fbl_ 


cof 


A 

fcl. 


cof 


fdl„ 


cof 


fel c 


cof 


ffl c ) 



b 



A 

• • • 

Then the value of A, A, B, B F, F by integration 

such that 

* dA f d^A _ 0 dA 

A - at = i ^2 ' A - i at 



A block diagram to compute all of the variables in the 
set of equations is presented in Fig. 2. 

In the computer program that is used for simulation 
all six equations for six degrees of freedom are used, but 
are interested in less than six degrees of freedom. The 
same program can be used by setting the coupling terms of 
non-used equations equal to zero and one in terms of princi- 
pal diagonal, e.g. only three degrees of freedom are used 
in this study, surge, sway and yaw, then all coupling terms 
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Pig. 2 . 

TRANSFER FUNCTION BLOCK DIAGRAM SIX DEGREES OF FREEDOM (Eqn. A) 



19 




TRANSFER FUNCTION BLOCK DIAGRAM SIX DEGREES OF FREEDOM (Ecjn.R) 
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Pig. 2 

TRANSFER FUNCTION BLOCK DIAGRAM SIX DEGREES OF FREEDOM (Eqn.C) 
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TRANSFER FUNCTION BLOCK DIAGRAM SIX DEGREE OF FREEDOM (Eqri.P) 
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Fig. 2 



TRANSFER FUNCTION BLOCK DIAGRAM SIX DEGREES OF FREEDOM (Eqn.E) 
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TRANSFER FUNCTION BLOCK DIAGRAM SIX DEGREES OF FREEDOM (Eqrt.F) 
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t 



are set equal to zero and unused terms on the principal 



diagonal equal to one, for example : 



aa 


ba 


0 


0 


0 


fa 




• • 

A 




J 1 


ab 


bb 


0 


0 


0 


fb 




• • 

B 




I 2 


0 


0 


1 


0 


0 


0 




• • 

C 




T 3 


0 


0 


0 


1 


0 


0 




• • 

D 




Z 4 


0 


0 


0 


0 


1 


0 




• • 

E 




I 5 


af 


bf 


0 


0 


0 


ff 




«• 

F 




*6 



and the left side of the unused equations are set equal 
to zero. 

With this program non-linear terms can be added such 
as, rudder deflection of waves and wind, etc., which will 
be done by adding all of these whose sum is IF^ e.g. 

1F 1 = KA^x Dr+KA 2 x Ds+KA 3 x Db+NA 

where Dr, Ds and Db are rudder deflection, canard de- 
flection etc. NA is the sum of all non-linear 

terms that effect the surge equation (X equation) . 

The program that will be used for solving these equa- 
tions is the "Continuous Systems Modeling Program" (CSMP) 
[Ref. 3] in which all constants are declared in the first 
section and then set the value of matrix for aa, ab, ac 

and so on (in program AAA is used for aa, AAB for 

ab ....... AFF for ff ) . In the initial section values of 

the COFACTORS aa, ba. . . .are determined. All of the CO- 
FACTORS and the subprogram VALUE is used to compute. This 
subprogram finds the value of the determinant of the 
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matrix. For all of the COFACTOR terms the element is set 
equal to one and the rest of the elements in that row and 
column are set equal to zero. For example, to find the 
value of COF.aa the following array is obtained: 



1 


0 


0 


0 


0 


0 


0 


bb 


cb 


db 


eb 


fb 


0 


be 


cc 


dc 


ec 


f c 


0 


bd 


cd 


dd 


ed 


fd 


0 


be 


ce 


de 


ee 


fe 


0 


bf 


cf 


df 


ef 


ff 



(In the computer program BAA is used for a^a^, GAA for 

a 2 a 2 **’)’ After the value of A and all cofactors are 

determined, the dynamic section is used to determine BAA, 

BAB GAA, GAB (if those terms contain variables) . 

In the dynamic section all variables that are functions 

of time are determined. The defining relations of the 

variables are also included in the dynamics section, i.e. 

• •• 

UDOT = ADDOT (U=A) , U=ADOT etc. XH,YH,ZH are deter- 

mined and are the vector terms, X,Y,Z whose origin is fixed 
on the earth (relative to the earth) . 
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